Introduction

I began creating a figure by taking the fGSEA results and coloring them by pathway types, e.g. glucose vs lipid metabolism, etc. in Illustrator. Preliminary results look really good, so I want to generate these graphs through R here.
I’ll run DESeq2 → fGSEA on all timepoints together, then every individual timepoints, then take Hallmark and KEGG data for WT SPF vs GF data to generate figures. If this works, I can do the same thing for other comparisons as well, e.g. SPF WT vs KO.

For this run, I grouped pathways into categories, some broad and others specific:

  • Metabolism
    • Oxidative phosphorylation
      • This was the original pathway that spurred this analysis that I noticed to be specifically upregulated in the GF condition all the time. I have some thoughts on why this may be that I comment on in conclusions. Oxidative phosphorylation refers to energy production that culminates in the electron transport chain.
    • Carbohydrate related pathways
      • Pathways that deal with carbohydrate substrate metabolism, such as gluconeogenesis, or pentose phosphate pathways. This doesn’t include things like TCA cycle since that is not carb-specific.
    • Glucose related pathways
      • Specifically glycolysis, gluconeogenesis, glycogenolysis, glycogenesis, or other glucose-specific pathways.
    • Lipid related pathways
      • This includes fatty acid, steroid, and bile acid metabolism.
    • Bile acid pathways
      • Pathways specific to bile acid metabolism.
    • Peroxisome related pathways
      • Pathways involving the peroxisome or that are specific to the peroxisome, e.g. PPAR.
    • Protein/AA related pathways
      • Metabolism of proteins, amino acids, or amino acid derivatives.
    • All metabolism pathways
      • This is literally just a combination of all the above pathways as well as miscellaneous metabolic processes, such as xenobiotic or vitamin metabolism.
  • Immune related pathways
    • This includes everything from innate to humoral to generic immune-specific pathways (e.g. IL6-JAK-STAT).
  • Cancer related pathways
    • I noticed cancer-related pathways are enriched in the presence of microbes, which led me to make this category. A minor note, in KEGG and Hallmark cancer pathways are specifically delineated, but GO biological processes does not label pathways by relation to disease. Instead, I broadly take tags involving differentiation, morphology, or cell growth. This topic brings another important point - there are databases specific to immune or cancer gene sets or even cellular location , and we might want to tap into those for deeper analysis in the future (for example, what cellular locations are most affected between microbial conditions?).

Since I’ve run this DESeq2 → fGSEA pipeline many times before in the same manner on this dataset, I will skip many of the visualizations/graphs for quality control I’ve generated in the past for this.

Setup

Libraries imported.

knitr::opts_chunk$set(echo = TRUE)
klippy::klippy(position = c("top", "right"))
reqpkg <- c("feather", "DESeq2", "dplyr", "magrittr", "genefilter", "AnnotationHub", "org.Mm.eg.db", "GO.db", "vsn", "pheatmap", "biomaRt", "curl", "RColorBrewer", "viridis", "fgsea", "tidyverse", "DT", "ggplot2", "ggpubr")
for (i in reqpkg) {
  print(i)
  print(packageVersion(i))
  suppressWarnings(suppressMessages(library(i, quietly=T, verbose=T, warn.conflicts=F,character.only=T)))
}
## [1] "feather"
## [1] '0.3.5'
## [1] "DESeq2"
## [1] '1.24.0'
## [1] "dplyr"
## [1] '0.8.5'
## [1] "magrittr"
## [1] '1.5'
## [1] "genefilter"
## [1] '1.66.0'
## [1] "AnnotationHub"
## [1] '2.16.1'
## [1] "org.Mm.eg.db"
## [1] '3.8.2'
## [1] "GO.db"
## [1] '3.8.2'
## [1] "vsn"
## [1] '3.52.0'
## [1] "pheatmap"
## [1] '1.0.12'
## [1] "biomaRt"
## [1] '2.40.5'
## [1] "curl"
## [1] '4.3'
## [1] "RColorBrewer"
## [1] '1.1.2'
## [1] "viridis"
## [1] '0.5.1'
## [1] "fgsea"
## [1] '1.10.1'
## [1] "tidyverse"
## [1] '1.3.0'
## [1] "DT"
## [1] '0.13'
## [1] "ggplot2"
## [1] '3.3.0'
## [1] "ggpubr"
## [1] '0.2.5'
theme_set(theme_pubr())
select <- dplyr::select
mart <- useDataset("mmusculus_gene_ensembl", mart=useMart("ensembl"))
mouseProteinCodingGenes <- getBM(attributes=c("ensembl_gene_id","external_gene_name","description"), filters='biotype', values=c('protein_coding'), mart=mart)
bm <- getBM(attributes=c("ensembl_gene_id", "hsapiens_homolog_associated_gene_name"), mart=mart) %>%
  distinct() %>%
  as_tibble() %>%
  na_if("") %>%
  na.omit()

pathway_folder <- '../pathways/'
f_res <- list()
anno <- list()
plot_res <- list()

dict <- list()
dict$immune <- c(
  "IMMUNE|INFLAMMATORY|TNFA|ALLOGRAFT|JAK_STAT|TGF_BETA|IL2_STAT5|INTERFERON_GAMMA|REACTIVE_OXYGEN_SPECIES|CYTOKINE|CHEMOKINE_SIGNALING|FC_EPSILON_RI|HEMATOPOIETIC_CELL|IMMUNODEFICIENCY|GRAFT_VERSUS_HOST|NOD_LIKE_RECEPTOR|TOLL_LIKE_RECEPTOR|FC_GAMMA_R|LEUKOCYTE|NATURAL_KILLER_CELL|T_CELL|COAGULATION|JNK_CASCADE|INTERLEUKIN|B_CELL|CELLULAR_DEFENSE_RESPONSE|LYMPHOCYTE|JUN_KINASE|CHEMOTAXIS|MACROPHAGE|MAST_CELL|MEGAKARYOCYTE|LYMPHOID_PROGENITOR_CELL|EXTRAVASATION|EOSINOPHIL|DENDRITIC_CELL|NEUTROPHIL|TUMOR_NECROSIS_FACTOR|RESPONSE_TO_BACTERIUM|DEFENSE_RESPONSE|COMPLEMENT_ACTIVATION|ANTIGEN|INTERFERON_ALPHA|T_HELPER|PLASMINOGEN"
  )
dict$cancer <- c(
  "CANCER|CARCINOMA|LEUKEMIA|REGULATION_OF_CELL_MORPHOGENESIS|VASCULAR_ENDOTHELIAL_GROWTH_FACTOR|CELL_GROWTH|MORPHOGENESIS|DIFFERENTIATION|PROLIFERATION|HOMEOSTASIS_OF_NUMBER_OF_CELLS|CELL_FATE_COMMITMENT|TISSUE_HOMEOSTASIS|RESPONSE_TO_GROWTH_FACTOR|UV_RESPONSE_DN|MELANOMA|DNA_REPAIR|DNA_DAMAGE_RESPONSE|DNA_DOUBLE_STRAND_BREAK|REGULATION_OF_CELL_DEATH"
  )
dict$ox_phos <- c(
  "OXIDATIVE_PHOSPHORYLATION|MITOCHONDRIAL_PROTON_TRANSPORTING_ATP_SYNTHASE_COMPLEX_ASSEMBLY|ELECTRON_TRANSPORT|RESPIRATORY_CHAIN|ENERGY_COUPLED_PROTON_TRANSPORT|ATP_BIOSYNTHETIC|NUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|NUCLEOSIDE_BISPHOSPHATE_METABOLIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|CELLULAR_RESPIRATION|MITOCHONDRIAL_MEMBRANE_POTENTIAL|CITRATE_CYCLE|TCA_CYCLE|TRICARBOXYLIC_ACID|POSITIVE_REGULATION_OF_ATPASE"
  )
dict$glucose <- c(
  "GLUCOSE|GLYCOLYSIS|GLUCONEOGENESIS|GLYCOGENESIS|GLYCOGENOLYSIS"
  )
dict$carb <- c(
  "GLUCOSE|GLYCOLYSIS|GLUCONEOGENESIS|GLYOXYLATE_AND_DICARBOXYLATE_METABOLISM|PYRUVATE_METABOLISM|CARBOHYDRATE_DERIVATIVE_BIOSYNTHETIC_PROCESS|MONOSACCHARIDE_BIOSYNTHETIC_PROCESS|CARBOHYDRATE_DERIVATIVE_METABOLIC_PROCESS|GLYCOSYL_COMPOUND_BIOSYNTHETIC|GLYCOSIDE_METABOLIC|GLYCOSYL_COMPOUND_CATABOLIC|RIBOSE_PHOSPHATE_BIOSYNTHETIC|GLYCOGENESIS|GLYCOGENOLYSIS|PENTOSE_PHOSPHATE_PATHWAY|CARBOHYDRATE_METABOLIC"
  )
dict$lipid <- c(
  "FATTY_ACID|CHOLESTEROL_HOMEOSTASIS|BILE_ACID|BIOSYNTHESIS_OF_UNSATURATED_FATTY_ACIDS|TERPENOID_BACKBONE_BIOSYNTHESIS|STEROID_BIOSYNTHESIS|ADIPOGENESIS|PHOSPHOLIPID_METABOLIC_PROCESS|LIPID_METABOLIC|LIPID_BIOSYNTHETIC|PHOSPHATIDYLCHOLINE_BIOSYNTHETIC|LIPID_HOMEOSTASIS|STEROL_METABOLIC|ACYLGLYCEROL_METABOLIC|STEROL_HOMEOSTASIS|CHOLINE_CATABOLIC|ISOPRENOID_METABOLIC|STEROID_METABOLIC|LIPOPROTEIN_METABOLIC|STEROID_BIOSYNTHETIC|LIPID_CATABOLIC|TRIGLYCERIDE_METABOLIC|LONG_CHAIN_FATTY_ACYL_COA_BIOSYNTHETIC|LIPID_OXIDATION|FATTY_ACYL_COA_BIOSYNTHETIC|CHOLESTEROL_BIOSYNTHETIC|STEROL_BIOSYNTHETIC|STEROID_HORMONE_BIOSYNTHESIS|DIACYLGLYCEROL_METABOLIC|GLYCOSPHINGOLIPID_BIOSYNTHETIC|TERPENOID_METABOLIC|LIPID_STORAGE"
  )
dict$bile <- c(
  "BILE_ACID"
  )
dict$peroxisome <- c(
  "PEROXISOME|PPAR_SIGNALING_PATHWAY"
  )
dict$protein <- c(
  "LYSINE_DEGRADATION|SELENOAMINO_ACID_METABOLISM|GLUTATHIONE_METABOLISM|TYROSINE_METABOLISM|ARGININE_AND_PROLINE_METABOLISM|CYSTEINE_AND_METHIONINE_METABOLISM|TRYPTOPHAN_METABOLISM|GLYCINE_SERINE_AND_THREONINE_METABOLISM|ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM|BETA_ALANINE_METABOLISM|PROTEASOME|VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION|REGULATION_OF_PROTEIN_CATABOLIC_PROCESS|PROTEOLYSIS|PROTEIN_CATABOLIC|PEPTIDE_METABOLIC|GLUTAMINE_FAMILY_AMINO_ACID_METABOLIC_PROCESS|TYROSINE_CATABOLIC|AMINO_ACID_CATABOLIC|AMINE_CATABOLIC|AMINO_ACID_FAMILY_METABOLIC|AMINO_ACID_METABOLIC|LYSINE_CATABOLIC|LIPOPROTEIN_METABOLIC|AMINO_ACID_BIOSYNTHETIC|GLUTAMATE_CATABOLIC|HOMOCYSTEINE_METABOLIC|L_PHENYLALANINE|AMINE_METABOLIC|PEPTIDE_BIOSYNTHETIC|PHENYLALANINE_METABOLISM|CARNITINE_BIOSYNTHETIC|ARGININE_BIOSYNTHETIC|GLUTATHIONE_METABOLIC|GLUTATHIONE_DERIVATIVE_BIOSYNTHETIC|TRYPTOPHAN_METABOLIC|PROTEIN_OXIDATION"
  )
dict$nucleic <- c(
  "PYRIMIDINE_METABOLISM|RNA_METABOLIC|RNA_CATABOLIC|NUCLEOSIDE_CATABOLIC|PYRIMIDINE_CONTAINING_COMPOUND_METABOLIC|PYRIDINE_CONTAINING_COMPOUND_METABOLIC|TRNA_THREONYLCARBAMOYLADENOSINE_METABOLIC|NUCLEOBASE_BIOSYNTHETIC|PYRIDINE_CONTAINING_COMPOUND_BIOSYNTHETIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_BIOSYNTHETIC|MRNA_CATABOLIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_METABOLIC|NUCLEOSIDE_PHOSPHATE_METABOLIC|PURINE_CONTAINING_COMPOUND_METABOLIC|NUCLEOSIDE_PHOSPHATE_BIOSYNTHETIC|PURINE_CONTAINING_COMPOUND_BIOSYNTHETIC|RIBONUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|CGMP_BIOSYNTHETIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_CATABOLIC"
  )
dict$scfa <- c(
  "BUTANOATE_METABOLISM|PROPANOATE_METABOLISM|BUTYRATE|ACETATE|GLUTAMATE|PROPRIONATE|ACETOACETIC_ACID"
  )
dict$metabolism <- c(
  "OXIDATIVE_PHOSPHORYLATION|GLUCOSE|GLUCOSE|FATTY_ACID|BILE_ACID|PEROXISOME|LYSINE_DEGRADATION|PYRIMIDINE_METABOLISM|BUTANOATE_METABOLISM|XENOBIOTIC_METABOLISM|MITOCHONDRIAL_PROTON_TRANSPORTING_ATP_SYNTHASE_COMPLEX_ASSEMBLY|GLYCOLYSIS|GLYCOLYSIS|CHOLESTEROL_HOMEOSTASIS|PPAR_SIGNALING_PATHWAY|SELENOAMINO_ACID_METABOLISM|RNA_METABOLIC|PROPANOATE_METABOLISM|LIMONENE_AND_PINENE_DEGRADATION|ELECTRON_TRANSPORT|GLUCONEOGENESIS|GLUCONEOGENESIS|BILE_ACID|GLUTATHIONE_METABOLISM|RNA_CATABOLIC|BUTYRATE|XENOBIOTIC_METABOLIC|RESPIRATORY_CHAIN|GLYCOGENESIS|GLYOXYLATE_AND_DICARBOXYLATE_METABOLISM|BIOSYNTHESIS_OF_UNSATURATED_FATTY_ACIDS|TYROSINE_METABOLISM|NUCLEOSIDE_CATABOLIC|ACETATE|METABOLIC|ENERGY_COUPLED_PROTON_TRANSPORT|GLYCOGENOLYSIS|PYRUVATE_METABOLISM|TERPENOID_BACKBONE_BIOSYNTHESIS|ARGININE_AND_PROLINE_METABOLISM|PYRIMIDINE_CONTAINING_COMPOUND_METABOLIC|GLUTAMATE|METABOLISM|ATP_BIOSYNTHETIC|CARBOHYDRATE_DERIVATIVE_BIOSYNTHETIC_PROCESS|STEROID_BIOSYNTHESIS|CYSTEINE_AND_METHIONINE_METABOLISM|PYRIDINE_CONTAINING_COMPOUND_METABOLIC|PROPRIONATE|CATABOLIC|NUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|MONOSACCHARIDE_BIOSYNTHETIC_PROCESS|ADIPOGENESIS|TRYPTOPHAN_METABOLISM|TRNA_THREONYLCARBAMOYLADENOSINE_METABOLIC|ACETOACETIC_ACID|CATABOLISM|NUCLEOSIDE_BISPHOSPHATE_METABOLIC|CARBOHYDRATE_DERIVATIVE_METABOLIC_PROCESS|PHOSPHOLIPID_METABOLIC_PROCESS|GLYCINE_SERINE_AND_THREONINE_METABOLISM|NUCLEOBASE_BIOSYNTHETIC|ANABOLIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|GLYCOSYL_COMPOUND_BIOSYNTHETIC|LIPID_METABOLIC|ALANINE_ASPARTATE_AND_GLUTAMATE_METABOLISM|PYRIDINE_CONTAINING_COMPOUND_BIOSYNTHETIC|ANABOLISM|CELLULAR_RESPIRATION|GLYCOSIDE_METABOLIC|LIPID_BIOSYNTHETIC|BETA_ALANINE_METABOLISM|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_BIOSYNTHETIC|BIOSYNTHETIC|MITOCHONDRIAL_MEMBRANE_POTENTIAL|GLYCOSYL_COMPOUND_CATABOLIC|PHOSPHATIDYLCHOLINE_BIOSYNTHETIC|PROTEASOME|MRNA_CATABOLIC|CELLULAR_RESPIRATION|CITRATE_CYCLE|RIBOSE_PHOSPHATE_BIOSYNTHETIC|LIPID_HOMEOSTASIS|VALINE_LEUCINE_AND_ISOLEUCINE_DEGRADATION|AEROBIC_RESPIRATION|TCA_CYCLE|GLYCOGENESIS|STEROL_METABOLIC|REGULATION_OF_PROTEIN_CATABOLIC_PROCESS|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_METABOLIC|RESPIRATORY_GASEOUS_EXCHANGE|TRICARBOXYLIC_ACID|GLYCOGENOLYSIS|ACYLGLYCEROL_METABOLIC|PROTEOLYSIS|NUCLEOSIDE_PHOSPHATE_METABOLIC|RESPONSE_TO_GLUCAGON|POSITIVE_REGULATION_OF_ATPASE|PENTOSE_PHOSPHATE_PATHWAY|STEROL_HOMEOSTASIS|PROTEIN_CATABOLIC|CARBOHYDRATE_METABOLIC|CHOLINE_CATABOLIC|PEPTIDE_METABOLIC|PURINE_CONTAINING_COMPOUND_METABOLIC|ISOPRENOID_METABOLIC|GLUTAMINE_FAMILY_AMINO_ACID_METABOLIC_PROCESS|NUCLEOSIDE_PHOSPHATE_BIOSYNTHETIC|STEROID_METABOLIC|TYROSINE_CATABOLIC|PURINE_CONTAINING_COMPOUND_BIOSYNTHETIC|LIPOPROTEIN_METABOLIC|AMINO_ACID_CATABOLIC|RIBONUCLEOSIDE_TRIPHOSPHATE_BIOSYNTHETIC|STEROID_BIOSYNTHETIC|AMINE_CATABOLIC|NUCLEOSIDE_BISPHOSPHATE_BIOSYNTHETIC|LIPID_CATABOLIC|AMINO_ACID_FAMILY_METABOLIC|CGMP_BIOSYNTHETIC|TRIGLYCERIDE_METABOLIC|AMINO_ACID_METABOLIC|NUCLEOBASE_CONTAINING_SMALL_MOLECULE_CATABOLIC|LONG_CHAIN_FATTY_ACYL_COA_BIOSYNTHETIC|LYSINE_CATABOLIC|LIPID_OXIDATION|LIPOPROTEIN_METABOLIC|FATTY_ACYL_COA_BIOSYNTHETIC|AMINO_ACID_BIOSYNTHETIC|CHOLESTEROL_BIOSYNTHETIC|GLUTAMATE_CATABOLIC|STEROL_BIOSYNTHETIC|HOMOCYSTEINE_METABOLIC|STEROID_HORMONE_BIOSYNTHESIS|L_PHENYLALANINE|DIACYLGLYCEROL_METABOLIC|AMINE_METABOLIC|GLYCOSPHINGOLIPID_BIOSYNTHETIC|PEPTIDE_BIOSYNTHETIC|TERPENOID_METABOLIC|PHENYLALANINE_METABOLISM|LIPID_STORAGE|CARNITINE_BIOSYNTHETIC|ARGININE_BIOSYNTHETIC|GLUTATHIONE_METABOLIC|GLUTATHIONE_DERIVATIVE_BIOSYNTHETIC|TRYPTOPHAN_METABOLIC|PROTEIN_OXIDATION"
  )
filter_fgsea <- function(selection, data, col1="#800000", col2="#D6D6CE") {
  on <- data %>% select(pathway) %>% filter(grepl(selection, pathway)) %>% pull()
  off <- data$pathway %>% setdiff(on)
  on <- data.frame(p=on, col=rep(col1, length(on)), in_category=rep(TRUE, length(on)))
  off <- data.frame(p=off, col=rep(col2, length(off)), in_category=rep(FALSE, length(off)))
  return(rbind(on, off))
}

create_palette <- function(x) setNames(as.character(x$p), x$col)

create_plot <- function(data, pal, title) {
  p <- ggplot(data, aes(x = reorder(pathway, NES), y = NES, fill = pathway)) +
    geom_col() +
    coord_flip() +
    labs(title = title) +
    theme(legend.position = "none", 
          plot.title = element_text(hjust = 0.5),
          axis.title = element_blank(), 
          axis.text.y = element_blank(), 
          axis.ticks.y = element_blank()) +
    scale_fill_manual(values = names(pal), limits = pal)
  return(p)
}

DESeq_run <- function(data) {
  condition_list <- data.frame(row.names=colnames(df), 
                               Genotype=rep(c(rep("WT",3), rep("KO",3)),2), 
                               Condition=c(rep("SPF",6),rep("GF",6)))
  condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
  condition_list$Condition <- relevel(condition_list$Condition, "SPF")
  dds <- DESeq2::DESeqDataSetFromMatrix(countData=data,
                                        colData = condition_list,
                                        design = ~ Genotype + Condition + Genotype:Condition)
  dds <- estimateSizeFactors(dds)
  dds$Group <- factor(paste0(dds$Genotype, dds$Condition))
  dds$Group <- relevel(dds$Group, "WTSPF")
  design(dds) <- ~  Group
  dds <- DESeq(dds)
  res <- results(dds,contrast=c("Group", "WTSPF", "WTGF"), tidy = TRUE)
  return(res)
}

All timepoints

This is analysis using raw featurecount data.

DESeq2

filePath = '../../../'
df <- read.csv(paste0(filePath, 'KF_RNASeq_counts_filtered.txt')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,] %>%
  set_rownames(.$Geneid) %>% select(-c(Geneid,X))

condition_list <- data.frame(row.names=colnames(df), 
                             Genotype=rep(c(rep("WT",18), rep("KO",18)),2), 
                             Condition=c(rep("SPF",36),rep("GF",36)))
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")

dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)
dds <- estimateSizeFactors(dds)
dds$Group <- factor(paste0(dds$Genotype, dds$Condition))
dds$Group <- relevel(dds$Group, "WTSPF")
design(dds) <- ~  Group
dds <- DESeq(dds)
WT_res <- results(dds,contrast=c("Group", "WTSPF", "WTGF"), tidy = TRUE)

fGSEA

res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
  select(hsapiens_homolog_associated_gene_name, stat) %>%
  na.omit() %>%
  distinct() %>%
  group_by(hsapiens_homolog_associated_gene_name) %>%
  summarize(stat=mean(stat))
ranks <- deframe(res)

f_res$time_all$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05) %>%
  arrange(desc(NES))

f_res$time_all$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))

f_res$time_all$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))

title <- ""

anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_all$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_all$hm <- lapply(pal, function(x) create_plot(f_res$time_all$hm, x, title))

anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_all$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_all$KEGG <- lapply(pal, function(x) create_plot(f_res$time_all$KEGG, x, title))

anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_all$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_all$GO <- lapply(pal, function(x) create_plot(f_res$time_all$GO, x, title))

Time 1 (ZT2)

DESeq2

filePath = '../../../data/R/by_time/'
df <- read_feather(paste0(filePath, 'time_1_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)

condition_list <- data.frame(row.names=colnames(df), 
                             Genotype=rep(c(rep("WT",3), rep("KO",3)),2), 
                             Condition=c(rep("SPF",6),rep("GF",6)))
condition_list$Genotype <- relevel(condition_list$Genotype, "WT")
condition_list$Condition <- relevel(condition_list$Condition, "SPF")

dds <- DESeq2::DESeqDataSetFromMatrix(countData=df,
                                      colData = condition_list,
                                      design = ~ Genotype + Condition + Genotype:Condition)
dds <- estimateSizeFactors(dds)
dds$Group <- factor(paste0(dds$Genotype, dds$Condition))
dds$Group <- relevel(dds$Group, "WTSPF")
design(dds) <- ~  Group
dds <- DESeq(dds)
WT_res <- results(dds,contrast=c("Group", "WTSPF", "WTGF"), tidy = TRUE)

fGSEA

res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
  select(hsapiens_homolog_associated_gene_name, stat) %>%
  na.omit() %>%
  distinct() %>%
  group_by(hsapiens_homolog_associated_gene_name) %>%
  summarize(stat=mean(stat))
ranks <- deframe(res)

f_res$time_1$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05) %>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.84% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_1$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.84% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_1$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.84% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 1"

anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_1$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_1$hm <- lapply(pal, function(x) create_plot(f_res$time_1$hm, x, title))

anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_1$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_1$KEGG <- lapply(pal, function(x) create_plot(f_res$time_1$KEGG, x, title))

anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_1$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_1$GO <- lapply(pal, function(x) create_plot(f_res$time_1$GO, x, title))

Time 2 (ZT6)

DESeq2

df <- read_feather(paste0(filePath, 'time_2_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)

fGSEA

res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
  select(hsapiens_homolog_associated_gene_name, stat) %>%
  na.omit() %>%
  distinct() %>%
  group_by(hsapiens_homolog_associated_gene_name) %>%
  summarize(stat=mean(stat))
ranks <- deframe(res)

f_res$time_2$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05) %>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.09% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_2$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.09% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_2$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.09% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 2"

anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_2$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_2$hm <- lapply(pal, function(x) create_plot(f_res$time_2$hm, x, title))

anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_2$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_2$KEGG <- lapply(pal, function(x) create_plot(f_res$time_2$KEGG, x, title))

anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_2$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_2$GO <- lapply(pal, function(x) create_plot(f_res$time_2$GO, x, title))

Time 3 (ZT10)

DESeq2

df <- read_feather(paste0(filePath, 'time_3_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)

fGSEA

res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
  select(hsapiens_homolog_associated_gene_name, stat) %>%
  na.omit() %>%
  distinct() %>%
  group_by(hsapiens_homolog_associated_gene_name) %>%
  summarize(stat=mean(stat))
ranks <- deframe(res)

f_res$time_3$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05) %>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_3$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_3$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 3"

anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_3$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_3$hm <- lapply(pal, function(x) create_plot(f_res$time_3$hm, x, title))

anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_3$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_3$KEGG <- lapply(pal, function(x) create_plot(f_res$time_3$KEGG, x, title))

anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_3$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_3$GO <- lapply(pal, function(x) create_plot(f_res$time_3$GO, x, title))

Time 4 (ZT14)

DESeq2

df <- read_feather(paste0(filePath, 'time_4_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)

fGSEA

res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
  select(hsapiens_homolog_associated_gene_name, stat) %>%
  na.omit() %>%
  distinct() %>%
  group_by(hsapiens_homolog_associated_gene_name) %>%
  summarize(stat=mean(stat))
ranks <- deframe(res)

f_res$time_4$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05) %>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.79% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_4$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.79% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_4$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.79% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 4"

anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_4$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_4$hm <- lapply(pal, function(x) create_plot(f_res$time_4$hm, x, title))

anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_4$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_4$KEGG <- lapply(pal, function(x) create_plot(f_res$time_4$KEGG, x, title))

anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_4$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_4$GO <- lapply(pal, function(x) create_plot(f_res$time_4$GO, x, title))

Time 5 (ZT18)

DESeq2

df <- read_feather(paste0(filePath, 'time_5_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)

fGSEA

res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
  select(hsapiens_homolog_associated_gene_name, stat) %>%
  na.omit() %>%
  distinct() %>%
  group_by(hsapiens_homolog_associated_gene_name) %>%
  summarize(stat=mean(stat))
ranks <- deframe(res)

f_res$time_5$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05) %>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.81% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_5$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.81% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_5$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (6.81% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 5"

anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_5$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_5$hm <- lapply(pal, function(x) create_plot(f_res$time_5$hm, x, title))

anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_5$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_5$KEGG <- lapply(pal, function(x) create_plot(f_res$time_5$KEGG, x, title))

anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_5$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_5$GO <- lapply(pal, function(x) create_plot(f_res$time_5$GO, x, title))

Time 6 (ZT22)

DESeq2

df <- read_feather(paste0(filePath, 'time_6_2.feather')) %>%
  .[rowSums(sapply(., '%in%', mouseProteinCodingGenes$ensembl_gene_id)) > 0,]
df <- set_rownames(df, df$Geneid) %>% dplyr::select(-Geneid)
WT_res <- DESeq_run(df)

fGSEA

res <- inner_join(WT_res, bm, by=c("row"="ensembl_gene_id")) %>%
  select(hsapiens_homolog_associated_gene_name, stat) %>%
  na.omit() %>%
  distinct() %>%
  group_by(hsapiens_homolog_associated_gene_name) %>%
  summarize(stat=mean(stat))
ranks <- deframe(res)

f_res$time_6$hm <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05) %>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "h.all.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_6$KEGG <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c2.cp.kegg.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
f_res$time_6$GO <- fgsea(pathways=gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), ranks, nperm=100000) %>%
  as_tibble() %>%
  filter(padj < 0.05)%>%
  arrange(desc(NES))
## Warning in fgsea(pathways = gmtPathways(paste0(pathway_folder, "c5.bp.v7.0.symbols.gmt")), : There are ties in the preranked stats (7.03% of the list).
## The order of those tied genes will be arbitrary, which may produce unexpected results.
title <- "time 6"

anno$hm <- lapply(dict, function(x) filter_fgsea(x, f_res$time_6$hm))
pal <- lapply(anno$hm, function(x) create_palette(x))
plot_res$time_6$hm <- lapply(pal, function(x) create_plot(f_res$time_6$hm, x, title))

anno$KEGG <- lapply(dict, function(x) filter_fgsea(x, f_res$time_6$KEGG))
pal <- lapply(anno$KEGG, function(x) create_palette(x))
plot_res$time_6$KEGG <- lapply(pal, function(x) create_plot(f_res$time_6$KEGG, x, title))

anno$GO <- lapply(dict, function(x) filter_fgsea(x, f_res$time_6$GO))
pal <- lapply(anno$GO, function(x) create_palette(x))
plot_res$time_6$GO <- lapply(pal, function(x) create_plot(f_res$time_6$GO, x, title))

Datatables

All timepoints

Hallmark

f_res$time_all$hm %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

KEGG

f_res$time_all$KEGG %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

GO biological processes

f_res$time_all$GO %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

Time 1

Hallmark

f_res$time_1$hm %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

KEGG

f_res$time_1$KEGG %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

GO biological processes

f_res$time_1$GO %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

Time 2

Hallmark

f_res$time_2$hm %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

KEGG

f_res$time_2$KEGG %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

GO

f_res$time_2$GO %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

Time 3

Hallmark

f_res$time_3$hm %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

KEGG

f_res$time_3$KEGG %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

GO

f_res$time_3$GO %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

Time 4

Hallmark

f_res$time_4$hm %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

KEGG

f_res$time_4$KEGG %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

GO

f_res$time_4$GO %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

Time 5

Hallmark

f_res$time_5$hm %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

KEGG

f_res$time_5$KEGG %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

GO

f_res$time_5$GO %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

Time 6

Hallmark

f_res$time_6$hm %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

KEGG

f_res$time_6$KEGG %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

GO

f_res$time_6$GO %>%
  select(-leadingEdge, -ES, -nMoreExtreme) %>%
  DT::datatable(filter = "bottom")

Graphs

All timepoints

Hallmark

ggarrange(plot_res$time_all$hm$immune, plot_res$time_all$hm$cancer, plot_res$time_all$hm$metabolism, plot_res$time_all$hm$ox_phos, plot_res$time_all$hm$glucose, plot_res$time_all$hm$carb, plot_res$time_all$hm$lipid, plot_res$time_all$hm$bile, plot_res$time_all$hm$peroxisome, plot_res$time_all$hm$protein, plot_res$time_all$hm$nucleic, plot_res$time_all$hm$scfa, ncol = 6, nrow = 2, 
          labels = c("immune", "cancer", "all metabolism", "oxidative phosphorylation", "glucose", "carb metabolism", "lipid metabolism", "bile", "peroxisome", "protein metabolism", "nucleic acid metabolism", "SCFA metabolism"), 
          font.label = list(size=10, face="plain"), hjust = 0.5, label.x = 0.5) %>% 
  annotate_figure(top = text_grob("Hallmark - all timepoints", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_all$KEGG$immune, plot_res$time_all$KEGG$cancer, plot_res$time_all$KEGG$metabolism, plot_res$time_all$KEGG$ox_phos, plot_res$time_all$KEGG$glucose, plot_res$time_all$KEGG$carb, plot_res$time_all$KEGG$lipid, plot_res$time_all$KEGG$bile, plot_res$time_all$KEGG$peroxisome, plot_res$time_all$KEGG$protein, plot_res$time_all$KEGG$nucleic, plot_res$time_all$KEGG$scfa, ncol = 6, nrow = 2, 
          labels = c("immune", "cancer", "all metabolism", "oxidative phosphorylation", "glucose", "carb metabolism", "lipid metabolism", "bile", "peroxisome", "protein metabolism", "nucleic acid metabolism", "SCFA metabolism"), 
          font.label = list(size=10, face="plain"), hjust = 0.5, label.x = 0.5) %>% 
  annotate_figure(top = text_grob("KEGG - all timepoints", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_all$GO$immune, plot_res$time_all$GO$cancer, plot_res$time_all$GO$metabolism, plot_res$time_all$GO$ox_phos, plot_res$time_all$GO$glucose, plot_res$time_all$GO$carb, plot_res$time_all$GO$lipid, plot_res$time_all$GO$bile, plot_res$time_all$GO$peroxisome, plot_res$time_all$GO$protein, plot_res$time_all$GO$nucleic, plot_res$time_all$GO$scfa, ncol = 6, nrow = 2, 
          labels = c("immune", "cancer", "all metabolism", "oxidative phosphorylation", "glucose", "carb metabolism", "lipid metabolism", "bile", "peroxisome", "protein metabolism", "nucleic acid metabolism", "SCFA metabolism"), 
          font.label = list(size=10, face="plain"), hjust = 0.5, label.x = 0.5) %>% 
  annotate_figure(top = text_grob("GO BP - all timepoints", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Every timepoint

Immune

Hallmark

ggarrange(plot_res$time_1$hm$immune, plot_res$time_2$hm$immune, plot_res$time_3$hm$immune, plot_res$time_4$hm$immune, plot_res$time_5$hm$immune, plot_res$time_6$hm$immune) %>% 
  annotate_figure(top = text_grob("Hallmark – immune", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$immune, plot_res$time_2$KEGG$immune, plot_res$time_3$KEGG$immune, plot_res$time_4$KEGG$immune, plot_res$time_5$KEGG$immune, plot_res$time_6$KEGG$immune) %>% 
  annotate_figure(top = text_grob("KEGG – immune", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$immune, plot_res$time_2$GO$immune, plot_res$time_3$GO$immune, plot_res$time_4$GO$immune, plot_res$time_5$GO$immune, plot_res$time_6$GO$immune) %>% 
  annotate_figure(top = text_grob("GO BP – immune", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Cancer

Hallmark

ggarrange(plot_res$time_1$hm$cancer, plot_res$time_2$hm$cancer, plot_res$time_3$hm$cancer, plot_res$time_4$hm$cancer, plot_res$time_5$hm$cancer, plot_res$time_6$hm$cancer) %>% 
  annotate_figure(top = text_grob("Hallmark – cancer", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$cancer, plot_res$time_2$KEGG$cancer, plot_res$time_3$KEGG$cancer, plot_res$time_4$KEGG$cancer, plot_res$time_5$KEGG$cancer, plot_res$time_6$KEGG$cancer) %>% 
  annotate_figure(top = text_grob("KEGG – cancer", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$cancer, plot_res$time_2$GO$cancer, plot_res$time_3$GO$cancer, plot_res$time_4$GO$cancer, plot_res$time_5$GO$cancer, plot_res$time_6$GO$cancer) %>% 
  annotate_figure(top = text_grob("GO BP – cancer", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

All metabolism

Hallmark

ggarrange(plot_res$time_1$hm$metabolism, plot_res$time_2$hm$metabolism, plot_res$time_3$hm$metabolism, plot_res$time_4$hm$metabolism, plot_res$time_5$hm$metabolism, plot_res$time_6$hm$metabolism) %>% 
  annotate_figure(top = text_grob("Hallmark – all metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$metabolism, plot_res$time_2$KEGG$metabolism, plot_res$time_3$KEGG$metabolism, plot_res$time_4$KEGG$metabolism, plot_res$time_5$KEGG$metabolism, plot_res$time_6$KEGG$metabolism) %>% 
  annotate_figure(top = text_grob("KEGG – all metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$metabolism, plot_res$time_2$GO$metabolism, plot_res$time_3$GO$metabolism, plot_res$time_4$GO$metabolism, plot_res$time_5$GO$metabolism, plot_res$time_6$GO$metabolism) %>% 
  annotate_figure(top = text_grob("GO BP – all metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Oxidative phosphorylation

Hallmark

ggarrange(plot_res$time_1$hm$ox_phos, plot_res$time_2$hm$ox_phos, plot_res$time_3$hm$ox_phos, plot_res$time_4$hm$ox_phos, plot_res$time_5$hm$ox_phos, plot_res$time_6$hm$ox_phos) %>% 
  annotate_figure(top = text_grob("Hallmark – oxidative phosphorylation", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$ox_phos, plot_res$time_2$KEGG$ox_phos, plot_res$time_3$KEGG$ox_phos, plot_res$time_4$KEGG$ox_phos, plot_res$time_5$KEGG$ox_phos, plot_res$time_6$KEGG$ox_phos) %>% 
  annotate_figure(top = text_grob("KEGG – oxidative phosphorylation", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$ox_phos, plot_res$time_2$GO$ox_phos, plot_res$time_3$GO$ox_phos, plot_res$time_4$GO$ox_phos, plot_res$time_5$GO$ox_phos, plot_res$time_6$GO$ox_phos) %>% 
  annotate_figure(top = text_grob("GO BP – oxidative phosphorylation", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Glucose

Hallmark

ggarrange(plot_res$time_1$hm$glucose, plot_res$time_2$hm$glucose, plot_res$time_3$hm$glucose, plot_res$time_4$hm$glucose, plot_res$time_5$hm$glucose, plot_res$time_6$hm$glucose) %>% 
  annotate_figure(top = text_grob("Hallmark – glucose metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$glucose, plot_res$time_2$KEGG$glucose, plot_res$time_3$KEGG$glucose, plot_res$time_4$KEGG$glucose, plot_res$time_5$KEGG$glucose, plot_res$time_6$KEGG$glucose) %>% 
  annotate_figure(top = text_grob("KEGG – glucose metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$glucose, plot_res$time_2$GO$glucose, plot_res$time_3$GO$glucose, plot_res$time_4$GO$glucose, plot_res$time_5$GO$glucose, plot_res$time_6$GO$glucose) %>% 
  annotate_figure(top = text_grob("GO BP – glucose metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Carbohydrate metabolism

Hallmark

ggarrange(plot_res$time_1$hm$carb, plot_res$time_2$hm$carb, plot_res$time_3$hm$carb, plot_res$time_4$hm$carb, plot_res$time_5$hm$carb, plot_res$time_6$hm$carb) %>% 
  annotate_figure(top = text_grob("Hallmark – carbohydrate metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$carb, plot_res$time_2$KEGG$carb, plot_res$time_3$KEGG$carb, plot_res$time_4$KEGG$carb, plot_res$time_5$KEGG$carb, plot_res$time_6$KEGG$carb) %>% 
  annotate_figure(top = text_grob("KEGG – carbohydrate metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$carb, plot_res$time_2$GO$carb, plot_res$time_3$GO$carb, plot_res$time_4$GO$carb, plot_res$time_5$GO$carb, plot_res$time_6$GO$carb) %>% 
  annotate_figure(top = text_grob("GO BP – carbohydrate metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Lipid metabolism

Hallmark

ggarrange(plot_res$time_1$hm$lipid, plot_res$time_2$hm$lipid, plot_res$time_3$hm$lipid, plot_res$time_4$hm$lipid, plot_res$time_5$hm$lipid, plot_res$time_6$hm$lipid) %>% 
  annotate_figure(top = text_grob("Hallmark – lipid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$lipid, plot_res$time_2$KEGG$lipid, plot_res$time_3$KEGG$lipid, plot_res$time_4$KEGG$lipid, plot_res$time_5$KEGG$lipid, plot_res$time_6$KEGG$lipid) %>% 
  annotate_figure(top = text_grob("KEGG – lipid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$lipid, plot_res$time_2$GO$lipid, plot_res$time_3$GO$lipid, plot_res$time_4$GO$lipid, plot_res$time_5$GO$lipid, plot_res$time_6$GO$lipid) %>% 
  annotate_figure(top = text_grob("GO BP – lipid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Bile

Hallmark

ggarrange(plot_res$time_1$hm$bile, plot_res$time_2$hm$bile, plot_res$time_3$hm$bile, plot_res$time_4$hm$bile, plot_res$time_5$hm$bile, plot_res$time_6$hm$bile) %>% 
  annotate_figure(top = text_grob("Hallmark – bile acid", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$bile, plot_res$time_2$KEGG$bile, plot_res$time_3$KEGG$bile, plot_res$time_4$KEGG$bile, plot_res$time_5$KEGG$bile, plot_res$time_6$KEGG$bile) %>% 
  annotate_figure(top = text_grob("KEGG – bile acid", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$bile, plot_res$time_2$GO$bile, plot_res$time_3$GO$bile, plot_res$time_4$GO$bile, plot_res$time_5$GO$bile, plot_res$time_6$GO$bile) %>% 
  annotate_figure(top = text_grob("GO BP – bile acid", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Peroxisome

Hallmark

ggarrange(plot_res$time_1$hm$peroxisome, plot_res$time_2$hm$peroxisome, plot_res$time_3$hm$peroxisome, plot_res$time_4$hm$peroxisome, plot_res$time_5$hm$peroxisome, plot_res$time_6$hm$peroxisome) %>% 
  annotate_figure(top = text_grob("Hallmark – peroxisome", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$peroxisome, plot_res$time_2$KEGG$peroxisome, plot_res$time_3$KEGG$peroxisome, plot_res$time_4$KEGG$peroxisome, plot_res$time_5$KEGG$peroxisome, plot_res$time_6$KEGG$peroxisome) %>% 
  annotate_figure(top = text_grob("KEGG – peroxisome", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$peroxisome, plot_res$time_2$GO$peroxisome, plot_res$time_3$GO$peroxisome, plot_res$time_4$GO$peroxisome, plot_res$time_5$GO$peroxisome, plot_res$time_6$GO$peroxisome) %>% 
  annotate_figure(top = text_grob("GO BP – peroxisome", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Protein metabolism

Hallmark

ggarrange(plot_res$time_1$hm$protein, plot_res$time_2$hm$protein, plot_res$time_3$hm$protein, plot_res$time_4$hm$protein, plot_res$time_5$hm$protein, plot_res$time_6$hm$protein) %>% 
  annotate_figure(top = text_grob("Hallmark – protein/amino acid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$protein, plot_res$time_2$KEGG$protein, plot_res$time_3$KEGG$protein, plot_res$time_4$KEGG$protein, plot_res$time_5$KEGG$protein, plot_res$time_6$KEGG$protein) %>% 
  annotate_figure(top = text_grob("KEGG – protein/amino acid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$protein, plot_res$time_2$GO$protein, plot_res$time_3$GO$protein, plot_res$time_4$GO$protein, plot_res$time_5$GO$protein, plot_res$time_6$GO$protein) %>% 
  annotate_figure(top = text_grob("GO BP – protein/amino acid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Nucleic acid metabolism

Hallmark

ggarrange(plot_res$time_1$hm$nucleic, plot_res$time_2$hm$nucleic, plot_res$time_3$hm$nucleic, plot_res$time_4$hm$nucleic, plot_res$time_5$hm$nucleic, plot_res$time_6$hm$nucleic) %>% 
  annotate_figure(top = text_grob("Hallmark – nucleic acid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$nucleic, plot_res$time_2$KEGG$nucleic, plot_res$time_3$KEGG$nucleic, plot_res$time_4$KEGG$nucleic, plot_res$time_5$KEGG$nucleic, plot_res$time_6$KEGG$nucleic) %>% 
  annotate_figure(top = text_grob("KEGG – nucleic acid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$nucleic, plot_res$time_2$GO$nucleic, plot_res$time_3$GO$nucleic, plot_res$time_4$GO$nucleic, plot_res$time_5$GO$nucleic, plot_res$time_6$GO$nucleic) %>% 
  annotate_figure(top = text_grob("GO BP – nucleic acid metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

SCFA metabolism

Hallmark

ggarrange(plot_res$time_1$hm$scfa, plot_res$time_2$hm$scfa, plot_res$time_3$hm$scfa, plot_res$time_4$hm$scfa, plot_res$time_5$hm$scfa, plot_res$time_6$hm$scfa) %>% 
  annotate_figure(top = text_grob("Hallmark – SCFA metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

KEGG

ggarrange(plot_res$time_1$KEGG$scfa, plot_res$time_2$KEGG$scfa, plot_res$time_3$KEGG$scfa, plot_res$time_4$KEGG$scfa, plot_res$time_5$KEGG$scfa, plot_res$time_6$KEGG$scfa) %>% 
  annotate_figure(top = text_grob("KEGG – SCFA metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

GO biological processes

ggarrange(plot_res$time_1$GO$scfa, plot_res$time_2$GO$scfa, plot_res$time_3$GO$scfa, plot_res$time_4$GO$scfa, plot_res$time_5$GO$scfa, plot_res$time_6$GO$scfa) %>% 
  annotate_figure(top = text_grob("GO BP – SCFA metabolism", face = "bold", size = 14),
                  left = text_grob("pathways", rot = 90),
                  bottom = text_grob("Normalized Enrichment Score"))

Methods

I generated the items in the categories by taking the lists of results and filtering key words and phrases specific to the category that can be used to pull all results containing said phrase. For example, for “glucose” one phrase was glycolysis.

To save time, after every analysis I checked only the unique results to sort into categories.

Conclusions

There is greater immune-related function in SPF mice which can be attributed to the essential role of microbes in immune development. Surprisingly, there is also increased cancer-related pathways in SPF mice as well, and GF mice have enrichment in tumor-preventing pathways, which suggests a role between microbes and oncogene expression. This may be related to the inflammation model of cancer, where any type of chronic immune response drives oncogene function, or this may highlight the relationship and balance between healthy immune functioning and oncogenesis.

Furthermore, there is broadly increased metabolic processes in GF mice accross all types of macromolecules. A minor note on the results, while it may look like protein metabolism is the most enriched for metabolism in GF mice, this may be more a reflection of the many amino acid related pathways that might all actually be centrally connected. Nonetheless, the extreme increase in host metabolism can reflect several possibilities: the host is compensating for reduced absorptive and metabolic capabilities in the gut due to the lack of microbes, or without microbes key regulatory components for metabolic homeostasis is dysregulated, and the mice are forced to rely on systemic metabolic increase to fulfill basic energy requirements, or core physiologic defects due to impaired developement in the GF condition, such as reduced mitochondrial function or brown fat is leading to upregulation of metabolic processes. With regard to the latter, in light of the decreased host immunity and cellular growth processes, an alternative hypothesis could be that the mice simply have more energy to spend in the absence of immune or heavy mitotic demands.

Digging deeper into the increased metabolic functioning of GF mice, it seems they rely more on lipid and protein sources for metabolism rather than glucose, or generally carbohydrates. Perhaps this reflects decreased carbohydrate digestion and absorption without bacteria, as other groups have shown bacteria are important for metabolite production from many starches and complex sugars.

In any case, oxidative phosphorylation specifically is the most enriched pathway in GF mice in nearly every case, thus, cellular respiration is, without a doubt, higher in GF mice. I believe this might be a culmination effect of increased metabolism accross the board — i.e. metabolism of all macromolecules are feeding into the electron tranport chain in the end to generate useable energy. It is important to note that these samples are liver samples, so I can’t comment on oxidative phosphorylation versus uncoupled electron transport energy generation in thermogenesis, but it is clear that mitochondria are playing a key role in regulating energy production of GF mice.

In the future, I would want to include another category for mitochondria, since many mitochondrial pathways showed up and are interesting, such as apoptotic or mitochondrial RNA expression pathways. Also, I want to take a look at the leading edge genes for important and recurring pathways, such as oxidative phosphorylation, and see what genes are most important in gene set enrichment. To confirm these results, metabolic cage data looking into energy usage and expenditure will be extremely insightful. Adding a high fat diet to this system would also be interesting to see how the carb-lipid-protein metabolic homeostasis might shift in the GF mice.

# ggarrange(plot_res$time_all$hm$immune, plot_res$time_all$hm$cancer, plot_res$time_all$hm$ox_phos) %>% 
#   annotate_figure(top = text_grob("Hallmark", face = "bold", size = 14),
#                   left = text_grob("pathways", rot = 90),
#                   bottom = text_grob("Normalized Enrichment Score"))

# f_res$time_1$hm %>% .[rowSums(sapply(., '%in%', setdiff(.$pathway, f_res$time_all$GO$pathway))) > 0,] %>%   select(-leadingEdge, -ES, -nMoreExtreme) %>%
#     DT::datatable(filter = "bottom")

Session info

## R version 3.6.3 (2020-02-29)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 18362)
## 
## Matrix products: default
## 
## locale:
## [1] LC_COLLATE=English_United States.1252 
## [2] LC_CTYPE=English_United States.1252   
## [3] LC_MONETARY=English_United States.1252
## [4] LC_NUMERIC=C                          
## [5] LC_TIME=English_United States.1252    
## 
## attached base packages:
## [1] parallel  stats4    stats     graphics  grDevices utils     datasets 
## [8] methods   base     
## 
## other attached packages:
##  [1] ggpubr_0.2.5                DT_0.13                    
##  [3] forcats_0.5.0               stringr_1.4.0              
##  [5] purrr_0.3.3                 readr_1.3.1                
##  [7] tidyr_1.0.2                 tibble_2.1.3               
##  [9] ggplot2_3.3.0               tidyverse_1.3.0            
## [11] fgsea_1.10.1                Rcpp_1.0.4.6               
## [13] viridis_0.5.1               viridisLite_0.3.0          
## [15] RColorBrewer_1.1-2          curl_4.3                   
## [17] biomaRt_2.40.5              pheatmap_1.0.12            
## [19] vsn_3.52.0                  GO.db_3.8.2                
## [21] org.Mm.eg.db_3.8.2          AnnotationDbi_1.46.1       
## [23] AnnotationHub_2.16.1        BiocFileCache_1.8.0        
## [25] dbplyr_1.4.2                genefilter_1.66.0          
## [27] magrittr_1.5                dplyr_0.8.5                
## [29] DESeq2_1.24.0               SummarizedExperiment_1.14.1
## [31] DelayedArray_0.10.0         BiocParallel_1.18.1        
## [33] matrixStats_0.56.0          Biobase_2.44.0             
## [35] GenomicRanges_1.36.1        GenomeInfoDb_1.20.0        
## [37] IRanges_2.18.3              S4Vectors_0.22.1           
## [39] BiocGenerics_0.30.0         feather_0.3.5              
## 
## loaded via a namespace (and not attached):
##   [1] snow_0.4-3                    readxl_1.3.1                 
##   [3] backports_1.1.5               Hmisc_4.4-0                  
##   [5] fastmatch_1.1-0               splines_3.6.3                
##   [7] crosstalk_1.1.0.1             digest_0.6.25                
##   [9] htmltools_0.4.0               fansi_0.4.1                  
##  [11] checkmate_2.0.0               memoise_1.1.0                
##  [13] cluster_2.1.0                 limma_3.40.6                 
##  [15] annotate_1.62.0               modelr_0.1.6                 
##  [17] prettyunits_1.1.1             jpeg_0.1-8.1                 
##  [19] colorspace_1.4-1              blob_1.2.1                   
##  [21] rvest_0.3.5                   rappdirs_0.3.1               
##  [23] haven_2.2.0                   xfun_0.12                    
##  [25] crayon_1.3.4                  RCurl_1.98-1.1               
##  [27] jsonlite_1.6.1                survival_3.1-11              
##  [29] glue_1.3.2                    gtable_0.3.0                 
##  [31] zlibbioc_1.30.0               XVector_0.24.0               
##  [33] scales_1.1.0                  DBI_1.1.0                    
##  [35] xtable_1.8-4                  progress_1.2.2               
##  [37] htmlTable_1.13.3              klippy_0.0.0.9500            
##  [39] foreign_0.8-75                bit_1.1-15.2                 
##  [41] preprocessCore_1.46.0         Formula_1.2-3                
##  [43] htmlwidgets_1.5.1             httr_1.4.1                   
##  [45] acepack_1.4.1                 farver_2.0.3                 
##  [47] pkgconfig_2.0.3               XML_3.99-0.3                 
##  [49] nnet_7.3-13                   locfit_1.5-9.2               
##  [51] labeling_0.3                  tidyselect_1.0.0             
##  [53] rlang_0.4.5                   later_1.0.0                  
##  [55] munsell_0.5.0                 cellranger_1.1.0             
##  [57] tools_3.6.3                   cli_2.0.2                    
##  [59] generics_0.0.2                RSQLite_2.2.0                
##  [61] broom_0.5.5                   evaluate_0.14                
##  [63] fastmap_1.0.1                 yaml_2.2.1                   
##  [65] knitr_1.28                    bit64_0.9-7                  
##  [67] fs_1.3.2                      nlme_3.1-144                 
##  [69] mime_0.9                      xml2_1.2.5                   
##  [71] compiler_3.6.3                rstudioapi_0.11              
##  [73] png_0.1-7                     interactiveDisplayBase_1.22.0
##  [75] ggsignif_0.6.0                affyio_1.54.0                
##  [77] reprex_0.3.0                  geneplotter_1.62.0           
##  [79] stringi_1.4.6                 lattice_0.20-40              
##  [81] Matrix_1.2-18                 vctrs_0.2.4                  
##  [83] pillar_1.4.3                  lifecycle_0.2.0              
##  [85] BiocManager_1.30.10           cowplot_1.0.0                
##  [87] data.table_1.12.8             bitops_1.0-6                 
##  [89] httpuv_1.5.2                  R6_2.4.1                     
##  [91] latticeExtra_0.6-29           affy_1.62.0                  
##  [93] promises_1.1.0                gridExtra_2.3                
##  [95] assertthat_0.2.1              withr_2.1.2                  
##  [97] GenomeInfoDbData_1.2.1        hms_0.5.3                    
##  [99] grid_3.6.3                    rpart_4.1-15                 
## [101] rmarkdown_2.1                 shiny_1.4.0.2                
## [103] lubridate_1.7.4               base64enc_0.1-3
 

By Sumeed Yoyo Manzoor

smanzoor@uchicago.edu